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Abstract 



The relation between the elastic wave equation for plane, isotropic bodies and an 
underlying classical ray dynamics is investigated. We study in particular the eigen- 
' frequencies of an elastic disc with free boundaries and their connection to periodic 

rays inside the circular domain. Even though the problem is separable, wave mixing 
between the shear and pressure component of the wave field at the boundary leads to 
an effective stochastic part in the ray dynamics. This introduces phenomena typically 
associated with classical chaos as for example an exponential increase in the number 
^ . of periodic orbits. Classically, the problem can be decomposed into an integrable 

part and a simple binary Markov process. Similarly, the wave equation can in the 
high frequency limit be mapped onto a quantum graph. Implications of this result 
' for the level statistics are discussed. Furthermore, a periodic trace formula is derived 

, from the scattering matrix based on the inside-outside duality between eigen-modes 

' and scattering solutions and periodic orbits are identified by Fourier transforming the 

spectral density. 
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^ . 1 Introduction 

In the beginning of the 20th century, Debye studied the density of vibrational modes in 
a solid body in the context of his work on the heat capacity. He found that the average 
density is in leading order proportional to the volume of the body times the third power 
of the frequency. Corrections to Debye's result were found later involving contributions 
due to the surface of the body [Q, |2). The density of eigenfrequencies of a solid body 
contains apart form these smooth terms also oscillatory contributions, which build up 
the discrete spectrum of individual eigenmodes. These oscillatory corrections have been 
studied intensively over the last decade or so in the context of the Helmholtz and the 
Schrodinger equation. In the high frequency limit they are known to be related to peri- 
odic orbits of an underlying classical dynamics, that is, the ray dynamics in a billiard in 
the former or the Hamiltonian dynamics of the corresponding classical system in the latter 
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case m, ^. It was in particular observed that different formulas apply when the classical 
dynamics is integrable Q or chaotic |^]. The relation between the wave equation and a 
related deterministic ray-dynamics is less obvious in elasticity. The wave equations are 
vectorial and different wave modes with differing wave velocities coexist. The notion of 
chaos or integrability needs to be reexamined here, which is the main purpose of this paper. 



In what follows we shall assume that the elastic body consists of an isotropic material. 
To reduce the dimensionality of the problem, we will furthermore consider only bodies of 
the form of a thin plate or an infinite rod with constant cross section. The vibrations in 
the plate or rod decouple then into two classes, the in-plane and the anti-plane vibrations 
1^]; (for plates, this is only true as long as the wavelength is much smaller than the thick- 
ness of the plate). The problem is thus reduced to two spatial dimensions. We will focus 
here on the in-plane vibrations for which the wave equation is still vectorial which makes 
it more complex than say the scalar Helmholtz equation. The wave field can be decom- 
posed into two polarizations, that is, pressure and shear waves, which have different wave 
speeds. The two polarizations couple at the boundary for physically relevant boundary 
conditions. An underlying ray dynamics emerging at high frequencies has similarly two 
types of rays traveling at different speeds and conversion between polarizations take place 
at the boundary. Ray conversion introduces a stochastic component into the dynamics 
and may lead to a large increase of possible ray trajectories compared to deterministic 
billiards for the same domain shapes 

We shall discuss mainly the case of a circular disc here, a separable problem due to 
the spherical symmetry. The case of elastic bodies without symmetries and fully chaotic 
classical ray dynamics has been discussed in |^ including a comparison with the semiclas- 
sical quantization of chaotic systems The scattering from two circular cavities in an 
elastic medium has been treated in The common idea is to write the spectral density 
as the trace of the Green's function which can in turn be expressed as sum over classical 
periodic orbits. We shall derive such a trace formula for the elastic disc and compare the 
results with the numerically calculated spectrum. We will furthermore show that the wave 
equation as well as the classical ray dynamics still posses a degree of 'randomness' due to 
the wave mixing at the boundary even though angular momentum is conserved. 

The quantum spectra of systems whose classical dynamics is chaotic has been found 
to follow random matrix theory originally developed in nuclear physics, see for example 
|ll|, H]. In the elastic case spectral statistics coinciding with RMT has been observed 
in experiments for rectangular and stadium shaped plates ||l3|, |l4|, Recently, spectra 
of graphs have also be shown to behave quite similar to chaotic systems |l^; they 
posses a trace formula for the spectral density and show random matrix statistics. We 
will make a connection between the ray dynamics in an elastic disc and the dynamics on 
a simple Markov graph, and will show how mode conversion effects the correlations in the 
eigenfrequency spectrum of the disc. 

The paper is organized as follows; we shall first introduce the elastic wave equation 
and a high frequency approximation of its boundary element kernel in section ^. Next, the 
classical ray dynamics in a disc is discussed in section |3[ In section ^, the exact solutions 
of the wave equation for circular symmetries is derived and high frequency approximations 
are discussed. We will then study the so-called nearest neighbor spacing distribution for 
the disc spectrum in more detail. An expression for the oscillatory part of the level density 
in terms of periodic orbits will be derived from the scattering matrix in section 0. 
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2 The elastic wave equation and short wavelength approxi- 
mations 



We shall consider the propagation of elastic deformations through an isotropic body. The 
partial differential equation in the frequency domain is the linear Navier-Cauchy equation 

fiA{u) + {X + fi)V(y-u)+puj^u = 0, (1) 

where u(x) is the displacement field in the body, A, /i are the material dependent Lame 
coefficients and p is the density. We shall restrict ourselves to two-dimensional problems in 
what follows. A generalization of the results in this section to three dimensions is, however, 
straightforward. The two-dimensional wave equation describes in-plane deformations in 
plates or wave propagation in cylindrical bodies extending to infinity along one axis. 

Introducing elastic potentials <I> and ^ by using standard Helmholtz decomposition of 
the displacement field u, that is, 

u = Up -I- Us with Up = V<i', Us = V X * , (2) 

the Navier-Cauchy reduces to two Helmholtz equations for the potentials 

(A + A:2)$ = 0; (A + A;2)* = 0. (3) 

Here, kp and kg are the wave numbers for the pressure (or longitudinal) and shear (or 
transversal) wave component, respectively. The wave velocities relating the wave numbers 
to the frequency iv via the dispersion relation kp^g = ^^/cp^s are different for the two different 
polarizations, on obtains 

j \ + 2p _ fjl 

Note that the pressure wave speed is always larger than the shear wave velocity. It is 
this difference in wave speed which leads to the phenomenon of wave splitting in the ray 
dynamics on impact with a boundary, see fig. |l[ The two wave equations (^) couple at 
the boundary, the details of the coupling depend on the boundary conditions. We shall in 
the following always consider free boundaries, that is, no forces act on the surface of the 
body. Forces acting on general surface elements are described in terms of the stress tensor 

aij = XdkUkSij + p {diUj + djUi) 

where the summation convention is used. Free boundary conditions correspond to 

t(u) = a{u) • n = (4) 

for the displacement field at the boundary where n denotes the normal to the boundary. 
The operator t refers to the traction. The traction operator (or traction matrix after 
representing it in a particular basis) for a circular boundary will be needed later to calcu- 
late the eigenfrequency spectrum of an elastic disc, and is explicitly derived in appendix^. 



Waves propagate freely inside the medium, that is, the pressure and shear component 
are decoupled and travel along straight lines. Wave splitting occurs at the boundary 
according to Snell's law 

Cp _ sin 9p 
c, sin^o' 
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Figure 1: Wave splitting for in-plane waves at the boundary 



where 9p, Op denote the angle of incident or reflection of the pressure and shear wave, 
respectively, measured with respect to the normal to the surface, see fig. ||. No mode 
conversion takes place for s-waves coming in at incident angles larger than a critical angle 
9c = arcsin (cs/cp). The reflection coefficients at impact with a plane interface for free 
boundary conditions can be given in terms of an orthogonal 2x2 coefficient matrix ot, 

sin 26s sin 26^ — cos^ 29s 
- - ^ (6) 



sin 29s sin 26*0 + k? cos^ 29, 



s 



ss 



Oips — otsp and a^pp ~t~ o-pg — 1 



where Otttt' relates an incoming wave of polarization vr G to an outgoing wave of 

polarization vr' and k = Cp/cg. In the literature, often only the reflection coefficients for 
the displacement field u are given P] related to the coefficient matrix a above by 



/ cos 
y Ctt' cos y^/ 

Here, lajr-n-'P is equivalent to the proportion of the energy density of the wave function 
undergoing transition from vr' — > vr whereas |a,r7r'P is the ratio of the corresponding energy 
fluxes normal to the boundary (with normal velocity CttCos^tt)- The unitarity of a thus 
implies flux conservation normal to the boundary. The tangential energy flux is, however, 
not conserved for free boundary conditions due to the non-vanishing tangential stress au 



giving rise to surface waves, (whereas Gnn = CTnt = crtn = at the boundary )||19|, gO|. 

We are interested here in solutions of the wave equations (||) in bounded domains in 
two dimensions. The set of eigenfrequencies is discrete and the solutions depend on the 
shape of the domain and on the boundary conditions. The wave equation is non-separable 
for typical domain shapes and numerical schemes, as for example boundary element meth- 
ods (BEM), have to be employed to calculate the eigenfrequencies and corresponding 
wave- functions. The BEM is typically less straightforward in the elastodynamical case 
compared to applying it to the scalar Helmholtz equation. The integral kernels become 
hyper-singular for common boundary conditions as for example free boundaries and the 
displacement field is vectorial. Standard techniques to apply BEM to the Navier-Cauchy 



equations are described in |21|. 



Relatively simple expression for the boundary integral kernels can, however, be ob- 
tained when considering the high frequency limit. A generalization of Bogomolny's transfer 
operator method p^], derived originally for the Helmholtz equation in bounded domains. 



4 



yields in two dimensions a boundary kernel in the form of a 2 x 2 matrix, that is, 



T{q,q';uj) 




sp 



fT^ikpL{q,q')-iUp^ Q 

fXgiksL{q,q')-iUs^ 



(8) 

where q, q' denote points on the boundary of the domain. L{q, q') is the distance between 
q and q' in the two-dimensional (x, y) - plane and the reflection coefficients are defined 
in @. The additional phases Up^s are Maslov indices which count the number of caustics 
along the path for each polarization Approximations to the eigenfrequencies are then 
obtained by solving 

det(l-T(w)) = 0. (9) 

The transfer operator (|8|) can be viewed as a discrete wave propagator acting on boundary 
wave functions by mapping outgoing two-component wave vectors at a point q on the 
boundary into outgoing wave vectors at q' . Wave mixing at the boundary enters through 
the matrix cx{q). Snell's law (^) is obtained by considering the two-step operator 



T\q, q"; uj) = j dq'T{q, q'; ij)T{q', q"; 



in stationary phase approximation. Considering n step operators T", one can derive 
periodic orbit trace formulas as presented in Q , see also section |5[ 

The transfer operator is in many respect a fairly crude approximation of the true BEM. 
It is the leading order term in a l/w expansions of the exact boundary integral kernel and 
does in particular not contain evanescent contributions. It can therefore not reproduce 
boundary effects like surface waves as well as diffraction or higher order mode mixing 
corrections. It is, however, a natural starting point to investigate the connection between 
the wave dynamics in elastic media of finite size and an underlying ray dynamics which 
includes ray splitting. 

In the following, we shall study the billiard with probably the most simple geometry, 
namely an elastic disc. Even though the wave equation is separable for this particular 
shape, there is some degree of wave chaos in this system, which can be traced back to 
chaotic components of an underlying ray dynamics. 



3 Classical ray dynamics for circular domains 

In this section, we shall discuss in more detail the ray dynamics in isotropic media of 
general shape and in particular for two-dimensional circular domains. 

We will adopt the following convention for a ray-trajectory in an elastic isotropic 
medium: a trajectory is at any instant in time given by its position and momentum as 
well as its polarization being either of p or s type. We may identify the wave numbers 
k-rr = cj/ctt, vt E {p,s} as the momenta of polarization vr with kg = nkp. A ray travels 
along straight lines between impacts with the boundary. At the boundary, it stays in a 
given polarization or undergoes mode conversion with probability 

ipp — ^ss — I'^^ppl ) ^ps — tsp — 1 (-'-O) 

where the reflection coefficients are given in (^). The angle of reflection depends on 
whether mode conversion takes place or not via Snell's law and so does the momentum. 
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A trajectory is uniquely determined after fixing its initial position and momentum and 
an infinite sequence of polarizations 7ri,7r2, . . ., vTj . . . G {s,p} reflecting the probabilistic 
nature of the dynamics. The dynamics in an elastic isotropic medium is thus taking place 
on two different energy sheets with energies Ep = kp and Eg = k1 with Ep < Eg. The 
dynamics on each sheet is deterministic, jumps from one sheet to the other may occur at 
the boundaries. 

We shall now turn to the dynamics in a circular disc with radius a. The angular 
momentum L is conserved at impact with a boundary both for rays staying in a given 
polarization and those undergoing mode conversion which follows directly from (|5|), that 
is, 

L = Zp sin 9p = Zg sin Og = const , (11) 

where we set Zp = akp and Zg = akg. The maximal angular momentum possible for fixed 
frequency uj is I-Lmaxl = Zg, no wave splitting occurs for Zp < \L\ < Zg. The dynamics in a 
disc follows a simple scaling relation and can be characterized in terms of the dimensionless 
impact parameters bj^ = L/zj^ with bp = nbg. The mode splitting regime is characterized 
by \bp\ < 1, pure s - rays exist for 1 < < k. A trajectory takes on at most two different 
angles of reflection 9p and 9g with sin^^r = 6,^- 

In the language of dynamical systems theory, one may say that the dynamics on each 
energy sheet is integrable, that is, trajectories for fixed L and a given polarization are 
confined to a two-dimensional manifold in phase space with the topology of a torus. Mode 
conversion couples two specific tori characterized by {Ep,L) and {Eg,L) and transitions 
are possible between these two tori only, see fig. ^ a. The total dynamics is thus not 
ergodic. 

For \bp\ < 1, the same initial condition in phase space does, however, lead to an 
exponentially increasing number of possible solutions due to the transitions between the 
energy sheets. This leads, for example, to an exponential increase in the number of 
periodic ray trajectories with increasing length, a phenomenon usually associated with 
classical chaos. The latter follows directly from the periodic orbit condition 

UpAipp + UgAifg = 27rm, Up + Ug = n > 2m . (12) 

Here, Aip.,^ is the change in the azimuthal angle for a ray with polarization tt between 
two reflections. The azimuthal angle and the angle of incidence are related through the 
relation Ayj^r = vr — 29-,^. The integer indices tIt^ correspond to the number of ray segments 
with polarization vr and m is the number of rotations around the center. The number 
N{n) of periodic orbits with n reflections (including permutations of the polarizations) 
increases thus exponentially like 

71 

N{n) ~ -2" . 

The fact, that these orbits (apart from the m = case) all form continuous families is, 
however, a feature known from integrable dynamics. 

The dynamics becomes particular simple when we restrict our considerations to the 
motion in the radial coordinate only. The dynamics in r is one-dimensional in each sheet 
and bound away from the center due to the centrifugal potential L^/2r^, see fig. |2| a. 
Transitions between the sheets occur at the boundary r = a for \bp\ < 1. The boundary 
map for the radial dynamics for fixed bp is thus a simple stochastic process which may be 
described in terms of a graph with two loops of the form shown in fig. |2| b. The transition 
rates (|lO|) which again depend only on the parameters bp and n define a Markov process 
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- boundary 



r = a r 



Figure 2: a) Radial dynamics in the disc takes place on two energy sheets 
Ep = kp and Eg = = K?Ep; transitions take place at the boundary r = a. 
b) The boundary map is equivalent to the probabilistic dynamics of a Markov 
process on a binary graph. 



on this graph with topological entropy hf = log 2 and exponential decay of correlation for 
ipp 7^ or 1. The chaotic component of the dynamics in the elastic disc is thus a two-level 
stochastic process. 

4 The elastic disc - exact results and the high-frequency 
limit 

The wave equation for a disc of radius a is separable independent of the boundary condi- 
tions. We will briefly discuss the exact solutions in the case of free boundaries and make 
a connection between the eigenfrequencies of the interior problem and the spectrum of 
the scattering matrix for the corresponding exterior problem. Details are referred to the 
appendices ^ and It will in particular be shown that the scattering matrix is equivalent 
to the transfer operator (^) in the high-frequency limit. 

4.1 The scattering matrix 

The elastic wave equation can for circular domains in two dimensions be solved in terms 
of the basis functions 

Mr,y^) = MKr)e''^, (13) 

where Ji{kT^r) is the lib. order Bessel function with vr = p or s and (p is the azimuthal angle. 
The separability of the wave equation reflects the conservation of angular momentum in 
the classical ray dynamics; we obtain as usual that the angular momentum takes on only 
integer values /. 

Applying free boundary conditions to a displacement wave vector obtained from the 
potentials (|l^) for fixed I leads to the secular equation, see appendix 0, 

det(tO(t^) = (14) 
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where the traction matrix t is given as 

f = f^^^ ~ h^s)Jlizp) - Zpj'iizp) il{Ji{zp) - Zpj'iizp)) \ 

' \ il{Jl{Zs) - Zsj[{Zs)) Zsj[{Zs)-{l^-\zl)Jl{Zs)) ^ ' 

and = kT^a. The condition (^) can be rewritten in terms of the scattering matrix for the 
outside problem, that is, for in-plane wave - scattering in an infinite plate with a circular 



hole. This connection is known as the inside-outside duality |23] between eigensolutions of 
the interior problem and transparent scattering solutions. The scattering matrix for fixed 
angular momentum / is given as [^] 

Si = -ti . (t+)-i , (16) 



see also appendix |B|, where t^, are obtained from the traction matrix (|15| ) by replacing 
the Bessel function and its derivatives by incoming and outgoing Hankel functions. Using 
Ji(z) = {hI^\z) + h['\z))/2 which implies the same identity for the corresponding 
traction matrices, that is, 

t/ = ^(tz+ + tr), (17) 



= det(tO = -det(t+)det(l + tr-(t+)-i) 



the eigenfrequency condition ([ij) can be written as 

1 

4 

= ^det(t+)det(l - S/). (18) 

The zeros of the first factor in ( [I^ ) are related to the resonances for exterior scattering 
at a circular cavity and are all in the lower complex z - plane. An eigenfrequency for the 
interior problem of the disc implies that the scattering matrix at the same frequency has 
an eigenvalue 1. That is, a scattering solution exists for which the obstacle, here the disc, 



is transparent. This principle holds for general shapes |23]. 



4.2 The scattering matrix in the high frequency hmit 

In the following we shall derive an approximation to the scattering matrix in the high 
frequency limit. 



The mode conversion regime |6p| < 1: In the energy - angular momentum regime 
for which the impact parameter = \L\/ Zp < 1, the Hankel functions entering may 
be split in terms of phases and amplitudes using the oscillatory Debye approximation; one 
obtains to leading order (appendix]^, 



with 




and cxi is the unitary matrix of reflection coefficients defined in (|6|) with angles of incidence 
fixed by the angular momentum condition ( pi] ) with L = I. We thus obtain for transitions 
between polarizations vr, tt' 

Si(7r^7r')~a^^/-e-^(<^-+'^-'). (21) 
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Note that the unitarity of the S - matrix is preserved in this approximation. This is 
not true in general for short wavelength approximations and results here from the quasi 
one-dimensionality of the dynamics. 



(22) 



Regime of no mode conversion 1 < < k : For angular momenta corresponding 
to incident angles larger than the critical angle or |6p| > 1, a somewhat different treatment 
needs to be employed. Here the exponential Debye expansion must be used for the pressure 
wave leading to an S - matrix of the form 

T 

with a reflection coefficient (in agreement with the plane interface result) 

o^ss = (23) 

and 

Zi = l + cos AOs + iS cos 9 s sin^ Os y^sin^ 9s - I/k^ . (24) 

Here, the boundary conditions lead to a pure phase shift dependent on the angle of inci- 
dence. There is no attenuation associated with this reflection contrary to the wave splitting 
case. The phase shift is due to a coupling to a surface longitudinal wave [18|. 



The transfer matrix for \bp\ < 1: By parametrising the boundary in terms of the 
azimuthal angle if, one obtains for the transfer operator (|8|) in a circular domain 

with A(y9 = \(p — (p'\ and d is the distance between two points 99, on the boundary, that 
is, 

d{Lp, ip') = 2a sin . 

The transfer operator depends only on the difference Aip and block-diagonalises with 
respect to the Fourier basis |/ >= ex.p{ilip) / ^/2tt for integer I. One obtains after evaluating 
the second integral by stationary phase 

<l\T\l'>=6u'ai Q ^2.<^J , (26) 

where the phases 

4>TT = ^Kd{A(p^) - l^Y' ~ 

are taken at the stationary phase point 

akT^ cos — ^ = a^Tr sin 9.,^ = I (28) 

which is the angular momentum condition (0). The phases in coincide with ( |20|) 
after inserting the stationary phase condition A(/?^/2 = arccos6,r, Eq. ([2^). The S-matrix 
in the approximation (|l^) is thus equivalent to the hermitian conjugate of the T-matrix 
up to a simple transformation in terms of a unitary diagonal matrix. In particular the 
eigenfrequency conditions det(l — T) ~ det(l — S) = coincide in the high frequency 
limit. 
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L 

Figure 3: The eigenspectrum of the disc for aluminum with k = 2.05; the 
wavenumber Zp = akp of an eigenfrequency is plotted versus the angular mo- 
mentum /. Eigenfrequencies with akp < I (or \bp\ > 1) are pure shear states, 
mode mixing occurs for akp > I. The lowest states in each / - series are due to 
surface waves (Rayleigh - waves). 



The mean density of eigenfrequencies for fixed /: The mean density of eigenfre- 
quencies di for fixed I can be obtained from the scattering matrix Si p3|, that is, 



Inserting the high frequency limit of the S matrix (|^) for \bp\ < 1 one obtains to leading 
order 

Id,, , , a 



diikp) = -— (0p + 0,) = - i^^l -hl + K^l - 62j . (29) 

Note that the mean level density depends on k-,^ and / only via the impact parameters h-,^. 
Equivalently, we obtain from (|2^) for 1 < |6p| < k 

Id, a 



IT dkp vr 



di{kp) = ——<Ps = -t^^l-hl. (30) 



4.3 Statistical properties of the eigenspectrum 

We saw in the last section that a transition takes place at \bp\ = \L\/kpa = 1 between 
a pure shear wave regime with \bp\ > 1 and a mode mixing regime with \bp\ < 1. This 
transition is reflected in the ray dynamics, see section Whereas only one family of 
(shear-) trajectories exists for fixed bp with \bp\ > 1, there are infinitely many such families 
for \bp\ < 1 and their number increases exponentially with the length of the trajectories. 
Such a phenomenon is reminiscent to the behavior typically found for chaotic classical 
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dynamics. It was observed in section |3| that the dynamics in the mode mixing regime can 
indeed be described by a stochastic Markov process on a two-loop graph, see fig. |2| b. In 
the hmit bp 0, however, the transition rates tpp = t^s approach one and the two modes 
decouple again leaving only two possible ray-trajectories. 

All these regimes should manifests itself also in the spectrum of the elastic disc. Spec- 
tral correlations are known to be particularly sensitive to the degree of chaos present in 
an underlying classical (ray-) dynamics ||Tl| , |l2|| . A popular statistical measure is the so- 
called nearest neighbour spacing distribution P{s) giving the probability of finding two 
adjacent eigenvalues of the spectrum (unfolded to mean level separation one) a distance 
s apart. P{s) follows a Poisson distribution for completely uncorrelated spectra, but has 
been conjectured to coincide with the results obtained for ensembles of random Hermitian 
matrices for completely chaotic dynamics. 

In fig. ^, the spectrum obtained from the exact eigenfrequency condition ( [l5| ) for an 
aluminum disc with k = 2.05 is shown. Here, the wavenumber Zp = akp of an eigenvalue 
with angular momentum / is plotted. One notices a difference in the mean density of 
eigenvalues for fixed / above and below the diagonal \hp\ = \l\/zp = 1, see Eqs. (|29|) 
and (|30[). The lowest eigenfrequency in each / - series can be attributed to a surface (or 
Rayleigh-) wave. Mode mixing occurs for \hp\ < 1 or akp > L. 

First we look at P{s) for the total spectrum which is obtained by projecting the dif- 
ferent I - eigenfrequency series in fig. ^ onto the Zp-axis. One finds indeed good agreement 
with Poisson statistics as shown in fig. ^ a. This reflects the fact that the ray-dynamics 
in the elastic disc is not ergodic but restricted to manifolds with fixed angular momentum 
L which have the form of a single torus for \bp\ > 1 or two coupled tori for \bp\ < 1. 
The eigenfrequencies for different / - series are thus uncorrelated, which leads to vanishing 
correlations in the full spectrum for large Zp after projection onto the Zp-axis. 



In order to see the influence of wave mixing on the spectrum, we need to study spectral 
correlations within a given / series. Indeed, the wave - dynamics for fixed I is given by the 
2x2 transfer matrix Eq. ( p^ ) (or the equivalent matrix for \bp\ > 1 obtained from (p2|)) 
being of the form 

The transfer matrix consisting of a unitary transition matrix a times a diagonal matrix 



is typical for propagation on quantum graphs as studied in |1(:, 17|. The phases in the 
diagonal matrix cpp, (pp are thereby interpreted as the lengths of the bonds in the graph, 
here the two loops in fig. ^ b. As we change Zp, we expect the correlations within a given 
/ series to change according to the degree of wave-mixing possible, that is, to the degree, 
that tpp deviates from zero or one. Note that the transition amplitude tpg allowing for 
transitions between p and s waves goes to zero both in the limit bp ^ and 6p — > 1, see Eq. 
(^; the wave modes decouple in these limits. We expect furthermore, that correlations 
are locally the same for different / series with bp fixed, that is, along straight lines in fig. Q 
with bp = const. We therefore study the nearest neighbour spacing distribution P{s) for 
eigenfrequencies lying in a window given by the intersection of the cone bp it Abp and the 
line L = I = const] the observed distributions P{s, bp) are indeed independent of / and we 
may average P{s,bp) over different / - series. 

The result is shown in fig. ^ b for three different values of bp and Abp = 0.01. In 
addition, the value of tpp as a function of bp is given; maximal mixing corresponds to 
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a) 



P(s) 




b) 



b = 0.12 









tpp 




Figure 4: Statistical properties of the eigenspectrum for aluminum with k = 2.05. 
a) the nearest neighbour spacing distribution (NNS) for all levels; b) NNS distri- 
butions for fixed impact parameter bp together with the transition probabilities 
tpp as a function of bp. The NNS is obtained for stretches of eigenvalues lying in 
a range bp it Abp for angular momenta from / = 300 — 3000; 



tpp = 0.5. The distributions are neither similar to the Poisson distribution, nor to any of 
the random matrix distributions. This non-universal behavior is typical for graphs with 
only a few bonds ||2^. One notices, however, that a gap is opening up in P{s) for small s 
as bp increases from zero. This can be interpreted as level-repulsion due to mode-mixing 
which increases as tpp deviates from one. For 6p — > 1 a different effect sets in; the pressure 
mode is suppressed and we witness the transition form a two mode to a one mode wave 
dynamics with equidistant eigenfrequencies. 

5 Eigenfrequency density and a periodic orbit trace formula 

So far we have shown that the traction matrix as well as the scattering and transfer matrix 
can be brought into block-diagonal form where each 2x2 block produces the spectrum for 
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fixed angular momentum /. In this section, we will make an explicit connection between the 
full spectrum and periodic trajectories in the elastic disc by looking at the total spectral 
density 

i 

where the sum runs over all eigenfrequencies (in terms of the angular velocity uj) of the 
disc. 

The spectral density can quite generally be written in terms of a smooth part and 
oscillatory contributions, the latter containing the period orbit contributions, that is, 

d{uj) = dsmoothi^) + dosc{^) ■ 

The smooth part gives the mean density of states which may be obtained using 

dN gYnooth 



dsmoothi^) 



where Ngmooth is the mean part of the spectral counting function N{uj) giving the number 
of levels below to. General results for the smooth part of the counting function of isotropic 
elastic media with free boundary conditions can be given I, g, that IS, 

Nsmoothi^j) ~ ^— — — Suj'^ + Luj + o{uj) (32) 

47r 47rCs 

with 

f3 = 4r]-3 + K + 4//7r, 

and 

' , P-'/^')^ 

1/k ' ... 



Here S, L are the surface area and the perimeter of the domain and finally 77 = cr/cs with 
cr the Rayleigh wave velocity 0. The leading term corresponds to the available phase 
space volume whereas the next order term contains corrections due to surface states. In 
contrast to the scalar Helmholtz equation for which the expansion of Ngmooth can be worked 
out to arbitrary order in principle |2^] only the first two terms are known in elastodynamics 
at present. 

The fluctuating part of the density of states which contains all the information about 
individual eigenfrequencies of the interior problem can be written in terms of the scattering 
matrix of the outside problem [23, 26]. One obtains 

dosciu^) = -Im ^ -— TV(S"(a;))t . (33) 
vr n duj 

n=l 

In the high frequency limit this term is related to periodic ray trajectories in the disc as 
will be shown below. 
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5.1 Oscillatory part of the density of eigenfrequencies 

We will first derive a periodic orbit expression for TrS" in the high frequency limit lo ^ 1 
using the block-diagonal form of the scattering matrix and the approximation (19), that 
is, 

TrS"(a;) = ^ TrSf 

K I ^^max 

« ^ ^^An e~^'^^'^p'^p~^""''^''K (34) 

Here, Imax ~ Zs denotes the maximal angular momentum and the second sum runs over 
all binary symbol strings 7r„ = 7ri7r2 . . . 7r„ of length n with VTj S {p, s}. The amplitude 
^7r„ is obtained as product over reflection coefficients (^), that is 

n 

1=1 

4>T^ are the phases defined in ( pO|) and rip, Ug equal the number of times the symbol p, s 
appears in tt^. The sum over binary symbol strings for fixed I is equivalent to a sum over 
all periodic paths in the binary graph fig. ||. Note that the reflection coefficients as well 
as the phases <&7r depend explicitly on / and lo. 

Next, we use Poisson summation to write the sum over I as 

TrS"' ~ A.-K e"^*'-"''"^'''''"''''^''-' 

TTn \l\<lmax 

OO plmax 

= / dM7r„e-2iK<^p+"»'^-)-2^^'^'. (35) 

7r„ m=-00 -^-^max 

Evaluating the integrals by stationary phase using 

= - arccos(//z,r) = , (36) 

where Ay^^r is the angle spanned by a ray-segment with polarization vr between two reflec- 
tions, we obtain the stationary phase condition 

^^totai = npAipp + UgAips = 27rm . (37) 



This is precisely the periodic orbit condition (ll2|), that is, only those angular momenta 
I* = cos ^Aip contribute significantly, for which a periodic orbit exists at frequency u>. 
The second derivative of the phases in (|3^) is 



_ djAiptotai - 27rm) ^ ^ / Up _^ Ug \ ^^^^ 
dl \ Zp cos 9p Zs cos 9s I 

with Ojr, the angle of incident and cos^tt = sin ^Aip-,^. 

After evaluating the phases in (35) at the stationary phase point I*, one obtains 
for the total phase 

TT vr 

2(np$p + Us^s) + 2-^1*171 = Upkpdp + Uskgds — n— = cuT — n— (39) 
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with n = Up + Us and d-j^ = 2a sin ^A^?, the length of a ray-segment of polarization vr 
between reflections. Furthermore, T is the period of the periodic orbit. We finally obtain 

(n) 
po 

where the sum is taken over all periodic rays with n reflections and 

Ari 



A, 



"^po 



+ 



By taking the complex conjugate and the derivative with respect to (jJ of eq. (|0|), we 
finally obtain the spectral density to leading order in 1/uj as 



I oo 1 (n) 

lauj X - 1 



d{Lj) w J — Y -^ApoTpo cos{ojTpo - mT/2 + tt/A) 



IT — ' n 

n=l po 




T A'' 
ppo _ ^ ^ cos{r{u^Tppo -nppo7r/2) + tt/A). (41) 



The last expression is obtained after summing over orbits related by cyclic permutations 
of the symbol code tt and the sum is now taken over all primitive periodic orbits (ppo) 
of arbitrary length, that is, over orbits not including repetitions and cyclic permutations. 
The second sum over r then includes the repetitions. 

We note in passing that the result (^Tj) can also be derived from a generalization of 
the abelian trace formula |27, S3] valid for systems with continuous symmetries. Here the 



symmetry is used to integrate over families of orbits F. Due to the rotational symmetry 
in the disc, one obtains 

/ 'Ji 71 \ 

^ orMtr^^vWWdL] V 2 47 

where Tr is the period of the orbit, a = fi — l with n the Maslov index, or the order of the 
(possibly) discrete symmetry group of the orbit and finally the product of the plane 
wave reflection coefficients for scattering at the boundary. Finally, 89 /dL, (also called the 
anholonomy matrix), describing the negative change of perimeter angle due to a change 
of impact parameter is 

— = +2 ^ + 



dL \ Zp cos dp Zs cos 9s 



in agreement with 
5.2 Periodic orbit spectrum 



Eqn. (41) gives an explicit connection between periodic ray trajectories in the disc and 
the eigenfrequencies of the system. By taking a Fourier transform of d{uj) one should be 
able to recover the periodic ray solutions including orbits which change polarization along 
their path. To suppress high-frequency oscillations in the signal, we convolute Eq. (41) on 
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both sides by a Gaussian test function as was also used in [^. The smoothing depends on 
a parameter r] proportional to the width, 

w{zp) = e-^^vh)\ Zp = kpa. (43) 

Figure ^ shows a comparison between a numerically calculated period spectrum ob- 
tained from the first 23000 eigenvalues of a disc using (|l^ and the approximative result 
([4l|), here for polyethylene with k = 3.61. The smooth part of the spectrum is removed 
and the density of eigenfrequencies is then Fourier transformed. This obtained period 
spectrum shows numerous peaks which fall roughly into three classes: orbits being of pure 
pressure, pure shear and mixed polarization type. In general the first class consists of 
the shortest orbits since the pressure waves have the fastest velocity. At f ~ 3.2 msec, 
we have an infinite number of pressure orbits accumulating at the boundary; orbits of 
higher winding number like the pentagram at t ~ 4.9 msec can also be resolved. Next, 
orbits with segments of both pressure and shear polarization type arise. Again one finds 
accumulation towards a limit orbit with the pressure segments becoming tangential to the 
boundary around t ~ 5.5 msec. We note here clear deviations of the actual numerical 
period spectrum from periodic orbit theory. Shear waves have incidence angles close to 
the critical angle and surface contributions become relevant. 

We note, that for the first time, periodic orbits changing polarization along their paths 
could clearly be identified in a Fourier spectrum of an eigenspectrum of an elastic body. 
Quite characteristic is the decay of peak height as the orbits increase in length. Note, 
however, that the class of pure shear orbits with long periods show up as comparatively 
high peaks in the spectrum. This is due to the lack of mode-conversion when the shear 
segments turn towards tangential incidence with a reflection coefficient becoming a pure 
phase (pSj). Surface orbits like the pure Rayleigh orbit at t « 12.2 msec can also be 
identified. 



6 Summary and Discussion 

We have studied the in-plane eigenfrequency spectrum for the elastic wave equation in 
two spatial dimensions with circular boundaries. It was shown that the eigenmodes can 
be expressed in terms of periodic rays of an underlying billiard-like classical dynamics. 
The ray dynamics conserves angular momentum, the wave equation becomes separable in 
polar coordinates, accordingly. It was pointed out, however, that the existence of two wave 
modes with different velocities partially destroys the integrability of the problem; the ray 
dynamics for fixed angular momentum takes place on two different energy manifolds in 
phase space for the two polarization. The dynamics on each manifold is one-dimensional 
and thus integrable, transitions between the energy sheets at the boundary introduce a 
purely probabilistic component. The classical dynamics corresponding to the elastic wave 
equation is therefore not deterministic and thus not integrable in the sense of Hamiltonian 
dynamics. 

By solving the wave equation explicitly and deriving high frequency approximations 
employing both the scattering matrix and the transfer operator, a connection between the 
wave dynamics in the disc for fixed angular momentum and the unitary propagation on 
a simple quantum graph could be established. Spectral correlations due to wave mixing 
manifest itself in a gap in the nearest neighbor spacing distribution and thus strong level 
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Figure 5: FFT of oscillating level density. The orbits are depicted with thin/fat 
lines for pressure/shear-polarization and are positioned according to their period. 
The periodic orbit theory ("p.o.t.") refers to (|4l|). The smoothing parameter in 
( |4^ is chosen here as r/ = 0.2. The actual material corresponds to polyethylene 
with Cp = 1950 m/s and Cs = 540 m/s. Finally the disc radius is a = Im. 
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repulsion. Finally, the full level density was expressed in terms of periodic orbits which 
could be identified explicitly in the Fourier transform of the exact density of eigenstates. 

The main corrections omitted in the high frequency approximations derived here occur 
for nearly tangential orbits both for pressure and shear components. This regime calls for 



a more refined approximation of the Bessel functions occurring in the traction matrix (15) 



as for example uniform approximations. Furthermore, periodic orbits accumulate at the 



boundary and the stationary phase approximation used to solve the integrals (35) breaks 
down. These corrections give rise to surface waves typical for free boundary conditions. 

Furthermore, higher order terms in the reflection coefficients become important when- 
ever the leading term Ops vanishes, that is, for normal incident 9 ^ 7r/2 corresponding to 
L — > 0. Periodic orbits at L = having both p and s segments can indeed been identified in 
the Fourier spectrum fig. ^. A detailed analysis of these effects will be discussed elsewhere. 
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A In-plane eigenfrequency spectrum for an elastic disc 

Due to the rotational symmetry of the disc, the wave equation separates in an angular 
and radial part. The eigenfunctions regular at the origin may be expressed in terms of 
Bessel functions; one writes the displacement field (§) as 

uj, = V(JKV)e"'''') ui = V X i{Ji{k,r)e~'''P) 

and p, s refers to pressure and shear polarization as usual, (also called "primary" and 
"secondary" wave in seismology referring to the time of arrival of these waves). A general 
interior eigenfunction can be expanded in these displacement fields, that is, 

= ai u|, + 02 . (44) 

To fulfill the free boundary condition, the traction of u/ has to vanish at the boundary, 
that is, 

t(ui) = ait(u|,)+a2t(ui) =0. (45) 
Writing the traction in terms of its radial and angular direction, one obtains 

= ait^(Up) + a2tr(u^) 

= ait<^(uj,) + a2t^(u^). (46) 
Expressed as a matrix equation this becomes 

(aia2) • = , (47) 
where we have collected the coordinates of both polarizations in a matrix 
tl = [t' ' 



7Vl. 

(/2 - ^Z^)Ji{Zp) -Zpj'iiZp) il{Jl{Zp) - Zpj[{Zp)) 

il{Jl{Zs) - Zsj'iiZs)) Zsj'liZs) - {t^ - \zl)Jl{Zs) 



(48) 
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with z = r or (/? and vr = p or s and thus ^7^1 = ^ii^ir)- We set as usual z-^ = ak-,^ with 
a, the radius of the disc. A superposition of these two polarizations fulfills the boundary 
condition only when 

det(ti) = (49) 

which is the eigenfrequency condition ([I^. For more details see for example [29| in which 
an expression equivalent to (Esl) is derived. 



B S- matrix 

B.l S-matrix in terms of traction operators 

Assume a general boundary condition 

n(u) = (50) 

with u an n-dimensional vector wave function and Q, a linear operator. Denote Uj the 
wave field being non-zero in its i-th component only with i € The scattering 

process of an incoming pure i-wave, u~, may then be described as 

3 

where denotes outgoing pure polarizations. The Eq. ( ^T[) must satisfy the boundary 
condition ( |50| ) for all i. Solving for the scattering matrix, one obtains 

S = -n(u~) • (n(u+))-i . (52) 

The operator ri(u) may be represented in matrix form 

{^{ui))j , (53) 

where the index i represents the component of the vector field u and j denotes the vector 
component of the operator ft, (that is, x, y or f, in two dimensions). 

In our case f2(u) is given by the traction t(u) (^) in polar coordinates. The analytic 
form is obtained from (^) using H\^\Hf'^ instead of Ji. 



B.2 The high frequency hmit ^ 1 for bp = l/zp < 1 

Starting from the Eq. (^) expressing the scattering matrix Si in terms of traction matrices, 
we will employ the oscillatory Debye approximation for the Hankel functions entering t^, 
that is, 

with ^ 

Q = yz^^—l? and = Q — 7r/4 — / arccos - . (55) 

z 

Here, z is either Zp or Zg with = a/cjr as usual. The traction matrices are obtained 
from ([l5|) by replacing the Bessel functions in terms of outgoing and incoming Hankel 
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functions, respectively. Inserting the Debye approximation and separating t in amplitude 
and phase, we obtain [^] 

Here 




= r i " 1 (57) 
and 

y+_{{l'-z'j2)-iQp nil -iQp) \ 

- [ a{l - tQs) iQs - {I' - zl/2) ) ' ^^^^ 

whereas Z~ has the same form as Z+ apart from replacing Q^^ by — Qtt- Hence, 

fe-'^" \ fe-'^" \ 



e"*'?''' y V e , 

with the unitary matrix cx defined as 

a = -G" • Z-/(G+ • Z+). (60) 

A straightforward but tedious expansion of a in terms of the wave numbers k^j^, expressing 
/ and in terms of the incidence angle 6-,^, that is, / = — sin 9.,^ and Q-^ = cos 9.^, one 
obtains indeed the reflection coefficients (^ in leading order. The limit taken corresponds 
to letting /ctt ^ oo and |/| — > oo but keeping their ratio fixed. As discussed above this 
ratio corresponds to fixed impact parameter /incidence angle. This finally reproduces 
formula (^). All the formulas are given here for negative angular momentum; choosing 
/ = ZTrSin^TT positive changes the sign of the off-diagonal components in (|l5|). This does 
not alter the eigenfrequency condition ( p!^ ) reflecting the degeneracy of the spectrum with 
respect to the sign change in /. 
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